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Two new formulations of general relativity are introduced. The first one is a parabolization of 
the Arnowitt, Deser, Misner (ADM) formulation and is derived by addition of combinations of 
the constraints and their derivatives to the right-hand-side of the ADM evolution equations. The 
desirable property of this modification is that it turns the surface of constraints into a local attractor 
because the constraint propagation equations become second-order parabolic independently of the 
gauge conditions employed. This system may be classified as mixed hyperbolic - second-order 
parabolic. The second formulation is a parabolization of the Kidder, Scheel, Teukolsky formulation 
and is a manifestly mixed strongly hyperbolic - second-order parabolic set of equations, bearing thus 
resemblance to the compressible Navier-Stokes equations. As a first test, a stability analysis of flat 
space is carried out and it is shown that the first modification exponentially damps and smoothes 
all constraint violating modes. These systems provide a new basis for constructing schemes for 
long-term and stable numerical integration of the Einstein field equations. 

PACS numbers: 04.25.Dm, 04.70.Bw 



I. INTRODUCTION 

For years people have tried to obtain analytic solutions 
of the complex field equations of Einstein's general theory 
of relativity (GR) . Apart from few cases where symmetry 
is invoked, it is almost impossible to analyze the compli- 
cated dynamics in the strong gravitational field regime as 
that is described by GR. Approximation methods have 
been developed over the course of time, but the most 
promising tool for tackling problems such as gravitational 
waves arising from binary black hole (BBH) or binary 
neutron star (BNS) mergers, gravitational collapse etc., 
is numerical relativity. 

In the past few years remarkable progress has been 
made towards achieving long term and stable evolution of 
the Einstein field equations. Recently particular cases of 
the longstanding problem of the evolution of BBH were 
solved [J, H, @. Despite this exceptional achievement, 
the general problem of long term and stable evolution of 
the GR equations remains open and there is still much 
work left to be done. There is no theory or prescription 
to chose what formulation(s) and what gauge conditions 
are suitable for the numerical solution of a given problem. 
For example, there does not seem to exist a definitive 
explanation of why the approaches of [J, H, 0] perform so 
well when contrasted to previous efforts. Furthermore, 
there are astrophysical and theoretical problems of great 
interest for which the formulations above have not been 
applied yet and it is not known whether they will prove 
successful in such cases. Such problems are the study 
of the internal structure of black holes and astrophysical 
phenomena, where except for black holes matter is also 
involved. 

The numerical integration of the Einstein equations is 
not an easy task, because the computations can become 
unstable and an exponential blow up of the numerical 
error may occur, even when the formulation employed 
admits a well-posed initial value problem. If the numer- 



ical techniques, gauge and boundary conditions chosen 
do not suffer by pathologies, perhaps the most impor- 
tant cause of instabilities during a free evolution is the 
growth of the constraint violating modes which are ex- 
cited by numerical errors. 

Several methods have been proposed to deal with this 
last type of instability. A path followed by many rela- 
tivists is referred to as constrained evolution [{J EU EH, 
EH, 0, EH EH , in which a subset of the dynamical vari- 
ables is evolved in time by use of the evolution equa- 
tions and the remaining variables are obtained by solv- 
ing the constraints after each time-step. This method 
has proven to be very accurate and has been successfully 
applied in cases where some type of symmetry is present. 
Another method, which resembles the approach of con- 
strained evolution, is called constraint projection [20ll2lj|. 
In this method all variables are evolved in time using 
the evolution equations and they are periodically "pro- 
jected" via an optimal algorithm to the closest configu- 
ration of dynamical fields which satisfies the constraints. 
This method is quite robust, but it can be computation- 
ally demanding because it requires the repeated solution 
of non-linear elliptic equations. 

Another strategy takes advantage of the fact that in 
the ideal case where the constraint equations are sat- 
isfied, one has the freedom to add combinations of the 
constraints to the right-hand-side (RHS) of the evolu- 
tion equations. This approach changes the mathematical 
properties of the equations, but not the physics these 
equations describe. A number of strongly and symmet- 
ric hyperbolic evolution equations have been constructed 
using this method, see [161 ] for a review, which resulted 
in improved stability of numerical simulations of GR. 

By virtue of the freedom to add the constraint equa- 
tions to the RHS of the evolution equations and/or by 
extension of the number of dynamical fields, it is also 
possible to add terms that act as "constraint drivers" , 
i.e., terms which turn the constraint surface into a local 



2 



attractor. This technique is nowadays usually termed as 
"constraint damping" . 

The constraint damping strategy naturally splits in 
two different categories: a) Additional, non-physical vari- 
ables are appropriately introduced in the evolution sys- 
tem, so that their time evolution drives the physical vari- 
ables toward the surface of constraints [H, EH , and b) no 
additional (non-physical) variables are introduced in the 
system, but only combinations of the constraints, which 
act as constraint drivers, are added to the RHS of the 
evolution equations of the physical fields, e.g. Q. In this 
case the constraint dissipation is inherent to the evolution 
equations of the physical variables themselves. 

The goal of many formulations in numerical relativ- 
ity has been to incorporate such drivers in symmetric or 
strongly hyperbolic systems without changing the princi- 
pal part of the evolution equations [H, EH . The common 
property of all these formulations is that the damping 
rate of the short-wavelength constraint-violating modes 
is independent of the wavelength of these modes. As a re- 
sult, the constraint violations are not smoothed during a 
free evolution, unless numerical dissipation is artificially 
added to a code. The existence of short-wavelength com- 
ponents may lead to spurious oscillations and potentially 
terminate the computations [J]. 

For this reason, it may be advantageous to consider 
formulations of GR which possess wavelength-dependent 
constraint damping. One way to achieve this goal is to 
construct formulations of GR that under free evolution 
force all or some of the constraints to evolve according 
to parabolic equations. Parabolic equations damp short 
wavelengths extremely efficiently and for this reason they 
posses desirable smoothing properties [29| . 

The purpose of this paper is to present two new pa- 
rameterized 3+1 formulations which fall in category b) 
of the constraint damping strategy. These two formula- 
tions will be referred throughout this work as system 1 
and system 2. The difference between these two systems 
and most formulations of GR, which belong to the same 
category, is that constraint dissipation does not only en- 
ter through low order terms, but also through higher or- 
der terms since the order of the of the evolution equations 
is increased. In system 1 and system 2 all or some of the 
constraint violating modes have been parabolized. 

The desirable property of the well-posed system 1 is 
that the evolution equations of the constraints become 
second-order parabolic, independently of the gauge con- 
ditions and spacetime configuration, which in turn im- 
plies that the constraint surface becomes a local attrac- 
tor. The evolution equations of system 1 may be classi- 
fied as a set of mixed hyperbolic- second-order parabolic 
partial differential equations (PDEs). 

The well-posed system 2, is a parabolization of what 
is usually referred to as the Kidder, Scheel, Teukolsky 
(KST) formulation which was introduced in [2| and later 
extended in Q . A generalization of the gauge condition 
chosen in [2j is given and a certain combination of the 
terms of system 1 is added to the RHS of the KST formu- 



lation. Unlike system 1, system 2 forces only a subset of 
the constraint variables to evolve according to parabolic 
equation and the constraint propagation system is mixed 
strongly hyperbolic - second-order parabolic (MHSP). 

It is here noted that hyperbolic and parabolic drivers 
have been constructed for the MHD equations, by intro- 
ducing additional non-physical fields, whose time evolu- 
tion drive the physical fields on the constraint surface, 
see for example [22j and references therein. A similar ap- 
proach has also been applied in numerical relativity via a 
Hamiltonian constraint relaxation method [23j , which in- 
troduces a parabolic equation to approximately solve the 
Hamiltonian constraint. The parabolic equation is solved 
not in real (coordinate) time, but in some fiducial time 
until the Hamiltonian constraint has been relaxed. The 
difference between the suggested approach in this work 
and the aforementioned two efforts is that here there is no 
need for extra variables and that the constraint damping 
is inherent to the evolution equations of the physical dy- 
namical variables, which implies that the damping takes 
place in real time and not in some fiducial time. In addi- 
tion to that, in contrast to the approach of [23|, system 
1 dissipates the perturbations which violate not only the 
Hamiltonian constraint, but all constraints of the formu- 
lation. 

Note also that a general, robust constraint driver has 
been developed and successfully applied to two formu- 
lations of Maxwell's equations in [2J|, where the evo- 
lution equations were cast into a type of MHSP set of 
equations. In [25], the same approach has been suc- 
cessfully applied to a "1+1" case of the ADM evolution 
equations, linearized about flat space, where the result- 
ing evolution system resembled the structure of a mixed 
hyperbolic - fourth-order parabolic set of PDEs. The 
approach of [24[ can in general be applied to the full 
non-linear equations of the first-order ADM formulation, 
which is re- introduced in this paper (overcoming this way 
the obstacle of fourth-order evolution systems). However, 
the main disadvantage of this method is that the result- 
ing evolution system would turn from quasi-linear into 
highly non- linear, see [25|, because the constraints would 
be added in quadrature. On the other hand, the ideas 
suggested in this work result in quasi-linear PDEs and 
this constitutes one of the major differences between this 
paper and the work in [24j . 

Finally, it must be mentioned that the idea of adding 
multiples of the derivatives of the constraints to the RHS 
of the evolution equations has been independently sug- 
gested in (2(|, where the approach was applied success- 
fully to different formulations of Maxwell's equations. 

This paper is organized as follows. In section HH a 
model equation is presented to illustrate the paraboliza- 
tion approach. The first order ADM formulation and 
gauges are presented in section HTT1 In section HVl system 
1 is presented and studied. A review of the KST formu- 
lation is presented, a generalization of the gauge of the 
original KST formulation is made and system 2 is given 
and studied in section [Vj In section I VII a study of the 
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stability properties of flat space with three different for- 
mulations is presented. In section I VIII a brief discussion 
of the choice of parameters is given. Finally, the work is 
concluded in section IVIIH 

II. A SIMPLE MODEL EQUATION 

To illustrate the strategy that will be employed in or- 
der to modify the Einstein equations a scalar wave equa- 
tion for the variable g(x,t) is introduced. A linear one- 
dimensional scalar wave equation is given by 

d t d t g = d x d x g, (1) 

where d x denotes a spatial derivative and dt a time 
derivative. Equation (TTJ) is equivalent to the following 
first-order system 

dtg =K, 

d t K =d x D, (2) 
d t D =d x K. 

Equations ([2]) have the same totality of solutions as equa- 
tion (fT]) only if the following constraint is satisfied at all 
times 

C = d x g - D = 0. (3) 

Equation ((3]) is essentially the definition of the D vari- 
able. In addition to the evolution equations of the dy- 
namical variables g, K, D, one can also obtain the evo- 
lution equation of the constraint C. By taking a time 
derivative of the constraint one finds 

d t C = 0. (4) 

The initial value problem in GR is similar to that of 
equations (J2|) and ((3|, in the sense that there is a set of 
evolution equations for the dynamical variables and a set 
of constraints that have to be satisfied at all times. In GR 
however, the equations are so much more complex that 
during their numerical integration the constraints may 
grow exponentially with time leading the computations 
to a termination. 

This is not the case with equation @, because © is 
a set of strongly hyperbolic PDEs (definitions and theo- 
rems pertaining to hyperbolic PDEs can be found in [29|) 
and hence a perfectly well-behaved one. It can be solved 
numerically, without any difficulty in satisfying the con- 
straint ([31) at all times. Therefore, the purpose of the 
model equation is to demonstrate a strategy in order to 
suppress very efficiently the growth of constraint violat- 
ing modes. 

This strategy is based on the fact that not only the 
constraint C has to be zero at all times, but also its spatial 
derivative. This implies that one can add combinations 
of those two to the RHS of the evolution equations 
without changing the solutions of the equations. With 



that in mind, the following modification is suggested in 
order to damp any constraint violating modes and still 
retain the well-poscdncss of the system. 

d t g =K + Xd x C, 

d t K =8 X D, (5) 
dtD ~d x K + (C. 

In order to prove that ([5]) admits a well-posed Cauchy 
problem, it is instructive to present from [29j some def- 
initions and theorems concerning parabolic and mixed 
parabolic-hyperbolic systems of quasi-linear partial dif- 
ferential equations . 

Consider a quasi-linear even-order operator of the 
form 

P(d x )=P 2m (d x ) + Q(d x ), (6) 

where 

P2m(d x )= ]T A v (u)d^d^---d d v \ (7) 

\u\—2m 

and 

Q(d x )= Mu)d^d 2 ^---d d ^, (8) 

\v\<2m-l 

where d v i = v i ^s" s denotes the partial derivative 

with respect to x a of order v s , s — 1, . . . , d with d the 
number of spatial dimensions, u is the column vector of 
the n dynamical variables and A u are n x n matrix func- 
tions of the dynamical variables, but not their deriva- 
tives. 

a. Definition 1. The equation dtu = P(d x )u 
is called parabolic, if for all neR d the eigen- 
values 0Jj(K,), j — 1,...,U of p2m(iK>) = 

_ Hv=2m A 2mKi Vl ■ ■ ■ n d Vd satisfy 

Re^in)) < -6\K\ 2m , j = l,...,n, (9) 
with some S > independent of k, where in the symbol 

b. Theorem 1 . The Cauchy problem for a parabolic 
system d t u — P(d x )u is well-posed. 

Assume now that a system of partial differential equa- 
tions can be written in the following form 

*(:)-(? (I «")(:)■ 

(. 10 ) 

where u and v are dynamical variables, P 2 = P2{d x ) is a 
second-order differential operator, P\ = P\{d x ) is a first- 
order operator and the coupling terms R\ 2 = R\2{d x ) 
and i?2i = R 2 i(d x ) are general first-order differential op- 
erators. Equations () 1 0[) may or may not have zeroth- 
order couplings. Further assume that the uncoupled sys- 
tems 

d t u = P 2 (d x )u 

d t v = P 1 (d x )v 1 J 

are second-order parabolic and first-order strongly hy- 
perbolic, respectively. 
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c. Definition 2. The coupled system (IT0|) is called 
mixed strongly hyperbolic - second-order parabolic. 

The coupling terms R±2 and R21 or any low order terms 
do not destroy the well-posedness of the system and the 
following theorem can be proved. 

d. Theorem 2. The Cauchy problem for mixed 
second- order parabolic- strongly hyperbolic systems (|10|) is 
well-posed. 

Equation {5} is a system of mixed strongly hyperbolic 
and second-order parabolic PDEs, provided that A > 0. 
This can be most easily seen, if the explicit form of ([5]) 
is written 

d t g = K + \d x d x g - Xd x D, (12) 
8 t K = d x D, (13) 
d t D = d x K + C(d x g-D). (14) 

It can be readily realized that equations (fT2|) -(fT4 | have 
the form of (fT0|) with u = g and v = {K, D}. The cou- 
pling terms are R12 = —Xd x D and R21 — (,d x g which are 
first order. Without these first-order couplings, equation 
([T2")) is second-order parabolic when A > 0, whereas the 
system of ([13")) and (|14p is strongly hyperbolic. There- 
fore, the total set of equations is MHSP and admits a 
well-posed Cauchy problem. 

The most interesting aspect of the modification intro- 
duced in iJSJ) is the structure of the evolution of the con- 
straint C. A straightforward calculation yields 

d t c = \d x d x c - cc. (15) 

Equation (fT5"]) is parabolic, provided that A > 0. There- 
fore, not only the evolution equations of the dynamical 
variables are well-posed, but also the evolution of the 
constraint. 

By Fourier transforming (|15|) , one can show that damp- 
ing of the constraint C takes place, if A > and ( > 0. 
This is because the damping rate of C is of the form 
exp(— (Ak 2 + C)t)i where k is the magnitude of the 
wavenumber. 

By the form of the damping rate of C it is evident that 
the value of £ can have a major impact on the long-term 
stability of the numerical integration of ([5]) , even though 
theorem 2 dictates that equation ([5]) admits a well-posed 
Cauchy problem irrespectively of the value of £. 

At this point it is instructive to try to understand why 
this is so. Note that the constraint and its derivative 
were added to the RHS of equation ([2]). Note also that 
the derivative of the constraint C may not be violated, 
while the constraint itself may be violated, as long as 
the violation occurs by a constant. Constant violations 
correspond to infinite wavelength and hence k = 0. This 
is the reason why there is no damping, if £ = and 
k = 0. To be more specific, at the RHS of the evolution of 
variable g the derivative of the constraint has been added. 
So, a violation of the constraint C by a constant is not 
"felt" by the A term. The A term damps only the finite 
wavelength constraint violating perturbations, because 
k =/= perturbations violate not only the constraint, but 



also the derivatives of the constraint. If £ < then in 
order for damping to take place, the following inequality 
has to be satisfied for all allowed n 

Xk 2 > -C (16) 

In a computer there are limits on the wavenumber of the 
perturbations, because of the finite length of the com- 
putational domain and the mesh size. The maximum 
k corresponds to the Nyquist frequency and the mini- 
mum to K m in = 27r/L, with L the length of the com- 
putational domain. Hence, condition (fT6|) is satisfied, if 
A > — C/K m in 2 - However, because of the efficiency of 
the damping very soon all allowed finite wavelength con- 
straint violating perturbations go to zero. As a result the 
A term in ([3J is switched off, since there are no further 
violations of the derivatives of the constraints. There- 
fore, the evolution of the constraint C is effectively given 
by equation (fl~5|) with A = 0. If in that case there is any 
residual k — violation of the constraints (which is in- 
evitable since an exact zero initial value of a constraint is 
never achieved in a computer), this violation will obtain 
explosive behavior (because £ < was assumed), which 
will occur at a time-scale Therefore, for long-term 
and stable numerical integration both A > and £ > 
are required. 

A final contradiction has to be resolved here. It was 
stated earlier that system ([5]) admits a well-posed Cauchy 
problem irrespectively of the value of £, as long as A > 0. 
However, in the previous paragraph it was argued that if 
£ < an exponential blow up of the numerical solution 
occurs. The resolution of this contradiction takes place, 
if one realizes that well-posedness does not in general 
guarantee global in time existence of solutions, but only 
short time existence. 

The difference between the suggested modification of 
this section and previous constraint damping strategies is 
the A term. The advantage of the new approach is that 
the shortest the wavelength of the constraint-violating 
perturbations the more efficiently those perturbations are 
getting damped. And it is the short wavelength perturba- 
tions, which produce the most spurious oscillations that 
may lead to blow up of the numerical error. Of course 
£ is important to damp the long wavelength (including 
k = 0) perturbations, which do not get damped very 
efficiently by the A modification. 

An important result of the wavelength-dependent 
damping timescale is that the constraint violating noise is 
smoothed with time. Another important property of the 
suggested modification is that it assigns non-zero "speed" 
of propagation to the static mode present i n (pt . This is 
may be an advantage because, according to [13], the pres- 
ence of zero eigenvalues can be responsible for numerical 
instabilities caused by perturbations by low order terms. 

The same parabolization strategy will be applied to 
the Einstein equations in section IIVI and it will be 
demonstrated that with appropriate modifications all 
constraints of the theory can be damped. 
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III. THE FIRST-ORDER ADM FORMULATION 
A. 3+1 decomposition of spacetime 

The four-dimensional spacetime of the theory of gen- 
eral relativity can be cast into a 3+1 decomposition 1]. 
This is achieved by assuming that the spacetime can be 
foliated by a one-parameter family of spacelike hypersur- 
faces. The four-dimensional metric is then written in the 
following form 

ds 2 = -a 2 dt 2 + ltj {dx l + P l dt){dx j + ftdt) (17) 

where jij is the positive definite 3-mctric on the t — 
const, hypersurfaces, a is the lapse function, (3 l the shift 
vector, and x % are the spatial coordinates, i = 1, 2, 3. 

To specify the spacetime foliation the gauge equations 
for the lapse function and the shift vector have to be 
given. These can generally be written as follows 

F a (^x b , a, P\d b a, d b d c a, ...,d b (3\ ...7^, cVfe, ...^ =0, 

a,b,c = 0,...,3, i,j = 1,2,3. 

_ (18) 
[28| . In this work only algebraic gauges of the following 
form are considered 



a = a(j), 



(19) 



where 7 is the determinant of the three- metric 7^, and 
the shift vector (3 % may be either constant or a fixed func- 
tion of the spacetime coordinates (i, x l ). In [27],[28| it was 
shown that (fT9|) makes the evolution well-posed on the 
surface of constraints, if 



A = dlna/dh^ > 0. 



(20) 



Working with (flT)|) , one has the flexibility to use either 
"1+log" slicing 



or a densitized lapse 



a = 1 + In 7 



(21) 



(22) 



where Q is constant or a fixed function of t, x % and a 
the densitization parameter. If (|2"Tj) is used, A = 1/a, 
whereas if (|22"|) is employed, A = a. Both gauge con- 
ditions have been used rather successfully in numerical 
relativity. However, note that "1+log" slicing has bet- 
ter singularity avoidance properties than the densitized 
lapse [301 ] . 



B. Evolution equations and constraints 

The Einstein equations in 3+1 split form are referred 
to as the ADM formulation and consist of two subsets 



of equations. The first subset is that of the evolution 
equations, which describe how the dynamical variables 
evolve in time. The second subset is that of the con- 
straint equations, which have to be satisfied for all times. 
The standard second-order ADM formulation [l| has as 
dynamical variables the 3-metric 7^ and the extrinsic 
curvature K%j of the 3D spacelike hypersurfaces. 

The simplest first-order form of the ADM formulation 
is derived from the second-order one by introducing ad- 
ditional dynamical variables 



kij 



dkiij, 



and then deriving the evolution equations for Dkij 
evolution equations of the first-order ADM are 



(23) 
The 



d Jij = -2aK, 



1 j > 



(24) 



d Kij = - ViVja + a {R l0 + KK l0 - 2 1 mn K im K ]n ) , 



d D k ij = —2adkKij - 2K i3 d k a, 



(25) 
(26) 



where K = j mn K mn is the trace of the extrinsic curva- 
ture, Vj is the covariant derivative operator associated 
with the three-metric, Rij is the Ricci tensor associated 
with the three-metric, and d — dt — £p, with £p the Lie 
derivative along the shift vector (3 l . 

The Lie derivatives of the dynamical variables are 



(27) 



£ l3 K ij = (ViP m )K mj + (V 3 l3 m )K mi + ^VmKij, 

(28) 



and 



£pDkij = (3 m d m Dkij + D mi jdk(3 m 

+ 2D km[l d j) p m + 2 lm(z d f) d k ^ 

The subset of constraint equations is 

H= R + K 2 - K mn K mn = 0, 



Mi = V m K r - 



V,A- = 0, i= 1,2,3, 



Ckij = dkjij — Dkij — 0. 



(29) 



(30) 



(31) 



(32) 



where R is the trace of the three-Ricci tensor. Ti is the 
Hamiltonian constraint, Mi the momentum constraints 
and Ckij are new constraints due to the introduction of 
Dkij- The equations as presented above do not include 
matter terms, because this work focuses on vacuum so- 
lutions of the Einstein equations. 

A final set of twelve constraints can be derived by 
virtue of (|23p . via taking a spatial derivative and keep- 
ing in mind that partial derivatives commute. Then Dkij 
must also satisfy 



■*klij 



d k D, 



diD. 



l-LJkij 



0. 



(33) 
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Constraints (|33| relate to (|32[) via the following equation 
Ckiij = diCkij — dkCuj. (34) 

It is important to note here that in equations (f2"4"l) - (f3"Tj) . 
all spatial partial derivatives of the three-metric are re- 
placed by the corresponding Dkij variable. However, in 
doing so an ambiguity arises due to the fact that there 
is not a single way to replace the mixed second-order 
derivatives of the three-metric by first order derivatives 
of the Dkij variables. For example the derivative did^riij 
can be replaced by either d\D2ij or c^-Diij. There is no 
rule or recipe on how to map those mixed derivatives of 
the standard ADM formulation to the first-order deriva- 
tives of the variables, when deriving the 1st order 
ADM formulation. Hence, one has to make a choice to 
lift this ambiguity. No matter what that choice is, in the 
continuum limit where Ckij = 0, there is no difference 
between the numerous final forms the first-order ADM 
formulation may obtain. However, in the discrete limit 
where the satisfaction of the Ckij constraints is at best 
approximate, there may be some difference. 

In this work the aforementioned ambiguity is lifted by 
choosing the most natural way in replacing the mixed 
derivatives of the three- metric, i.e., they are written in 
as most symmetric way as possible. This is done by fol- 
lowing the suggestion of Q , according to which 

dkdigij = d {k D l)ir (35) 

This last modification is equivalent to arbitrarily re- 
placing partial derivatives of 7^ by Dkij in the spatial 
Ricci tensor Rij , and then recalculating the Ricci tensor 
as follows 

Rij = Rij + -^l ls (Cujs +Cju s — Cijis), (36) 

where Rij stands for the recalculated Ricci tensor. It 
is noted here, that it seems that the last term in the 
parenthesis makes the Ricci tensor asymmetric under 
the interchange of its indices, because that term satisfies 
Cijki = —Cjiki- However, a straightforward calculation 
shows that the addition of this term is the only way that 
the Ricci tensor can be cast into a manifestly symmetric 
object, in the context of the 1st order ADM formulation. 

Note that a more general Ricci tensor R^ can be cal- 
culated by adding to the RHS of (|3"6")l extra terms as 
follows 

Rij =i?„+^ 7 Qb C afo)b . (37) 

This modification is redundant for system 1 because the 
system is well-posed and has constraint damping when 
if> = 0, see section ITVl However, ip 7^ is essential for the 
well-posedness of system 2, see section IVl 

C. Evolution of the constraints of lst-order ADM 

By taking the time-derivative, d Q , of the constraint 
equations of the first order ADM formulation one can 



derive the evolution equations for those constraints by 
replacing the time derivatives of the dynamical variables 
through the evolution equations and then eliminate the 
appropriate combinations of spatial derivatives of the dy- 
namical variables via the constraint equations. After 
some tedious algebra, one finds that the evolution of the 
constraints, without using (|36p. is given by 

d C ki j = 0, (38) 

d C mj = 0, (39) 

£) H =2aKH - 2 ai ks d k M s - AM nl nm d m a 
+ 2aY nn T s mn M s + H mns C mns , 

doMi = - ^adiH + ^a 7 mn Y s diCrmsn - ma 

+ aKM l + ^M msn C lmsn + { ^Mr snl C msnl 

, (3)i\rmnn. , (4) itrmsnp 

(41) 

where H mns K x )M msn , t 2 )^™"^ (3)j|f»» i (4)j^.msm 
are functions of the dynamical variables and their spatial 
derivatives, as well as of the lapse function and its spa- 
tial derivatives. Their exact forms are long and they are 
not presented here for brevity, but the full equations are 
available by the author upon request. 

The Lie derivatives of the constraints along the shift 
vector are 

£f}C M j = 0, (42) 

£pCkiij = P a d a Ckhj + 2(d[i(3 a )C k \aij + 2C k ia(j(dj)P a ), 

(43) 

£ p H = f3 m d m H, (44) 

£pM l = P m d m Mi + M m a,/3 m . (45) 

From equations ([38l) - (|4"Tj) it is expected that if the con- 
straints are satisfied on the initial time-slice, they should 
be satisfied throughout the entire evolution. 

It is noted here, that if equation ([56]) is used for the 
Ricci tensor, it is straightforward to show that the Ricci 
scalar and hence the Hamiltonian constraint remain un- 
affected. However, since the evolution equations of the 
extrinsic curvature are modified, the resulting expres- 
sions of the evolution of the momentum and Hamiltonian 
constraints change. Their exact expressions are not pre- 
sented here, but their derivation is straightforward, if one 
uses equations (|A1|) - (|A6|) of appendix |A1 

The principal part of (|38)) - (j4T|) . has only two non-zero 
eigenvalues, which propagate with the fundamental speed 
±a, when (3 l = 0. All remaining eigenvalues are zero 
(static modes). Static modes are present in the ADM 



7 



constraint evolution system, even if (3 % ^ because the 
time derivative of C k ij is zero independently of the shift 
vector. 

This is where the argument in jl7] fits in, to explain 
why several formulations of GR blow up when integrated 
numerically. Just like the ADM constraint propagation 
equations, the evolution of the constraints with other 
formulations, such as KST, have a set of eigenvalues. 
This implies that when constraint violating perturbations 
arise during numerical integration of the evolution equa- 
tions, these perturbations remain in the system because 
they cannot propagate off the computational domain. As 
a result, continuous perturbations by numerical errors 
lead to accumulation of constraint violating noise which 
eventually terminates the computations. 

Note, however, that the reason why the ADM formula- 
tion has not been successful in multi-dimensional simula- 
tions is because the evolution system is ill-posed in more 
than one dimensions. 

This section is concluded by noting that since the ADM 
constraint equations have to be satisfied at every single 
point of a spacelike hypersurface, not only the constraints 
themselves must be zero, but also their spatial partial 
derivatives of any order have to be zero, i.e., 



ds^Ckiij =0, 
d. v H =0, 
d s v Mi =0. 



(46) 



This means that one is free to add to the RHS of the 
ADM equations not only combinations of the constraints 
themselves, but also combinations of their spatial deriva- 
tives. This property, will be used in subsequent sections 
to obtain parabolized formulations of general relativity. 



IV. PARABOLIZATION OF THE 1ST-ORDER 
ADM FORMULATION 

In this section equations to (j2"6")l are modified in 
such a way as to have as many parabolic modes as pos- 
sible. Parabolic modes are modes which decay expo- 
nentially with time at a rate exp (— en 2 t), with e > 
a damping parameter and k 2 the magnitude squared of 
the wavevector. With that in mind, the following modi- 
fication is proposed 



dtHj = {ADM) + aX^dtCa 



(47) 



dtKij = {ADM) + (/xryijj^daMb + 9ad (l M 3 ) , (48) 



d t D kij = {ADM) + e ai ab d a C bktJ + ta 1%J d k H + ^C kij , 



(49) 



where {ADM) denotes the RHS of equations 
Equations (l47p ~ (H9^) are called system 1 in this work 



The added derivatives of the Hamiltonian and momen- 
tum constraints are given by rather complicated and long 
expressions and for this reason they are not presented 
here. However, their full expressions are available by the 
author upon request. It is only noted that in those ex- 
pressions whenever derivatives of the 3-metric, d k jij, are 
present they are replaced by the variables D k ij. This 
ways one ensures that the resulting system of PDEs is 
quasi-linear instead of non-linear. 



A. Structure and Well-Posedness of system 1 

Any new formulation of GR intended for numerical 
simulations is required to admit a well-posed initial value 
problem. This is because an ill-posed formulation cannot 
lead to long-term and stable numerical simulations of the 
Einstein Equations. This sub-section establishes the well- 
posedness of system 1 and discusses the structure of its 
evolution equations. 

In the analysis that follows C is considered to be zero 
because the results with ( ^ are rather complicated. 
It has been found that, just like in the model equation 
of section UH £ has no impact on the well-posedness of 
system 1, but it must satisfy £ > for long-term and 
stable numerical simulations. 

For the study of well-posedness of a given set of evo- 
lution equations the zeroth-order terms are omitted and 
only the higher-order terms have to be considered. For 
the family of gauges (|19[) the highest order operator of 
system 1 is of second-order and according to 29] it has 
to be analyzed first. 

The second-order part is given by the following equa- 
tions 



d t K i3 « a<f)-i lj -i mn -i ks {d m d k K ns - d k d s K mn )+ 



(50) 



-a9j mn {d m d z K jn + d m d 3 K in - 2d l d ] K mn ), 

(51) 



d t D Hj « a 7 £s 7 m "(a fe 9^ ms „ - d k d e D smn )+ 
aej mn {d m d n D klj ~ d m d k D nij ), 



(52) 



where the symbol "~" means equal to the second-order 
part. It can be readily seen that the evolution equations 
of 7-jj, Kij and D k ij decouple to second order. This im- 
plies that one can study each set of equations separately. 

The evolution equations for the three-metric (|50[) are 
quasi-linear second-order parabolic, if A > 0. All six 
eigenvalues of P{in) are equal to 

u>i = — a\n 2 , i = 1, ... 6 (53) 

and for A > condition ^ is satisfied. 
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The set of eigenvalues of P(in) of equation ([51]) is 

uji = — —a9n 2 i — 7, 8, 

ujg = 2acj)K 2 , V*Q 
Wi = 0, i = 10, . . . , 12. 

Therefore, -ftTjj does not evolve according to parabolic 
equations because condition ([9]) cannot be satisfied due 
to the existence of zero eigenvalues. 

An analogous analysis of equations ([52")) yields the fol- 
lowing eigenvalues 

Ui = — aen 2 , i = 13, . . . 24, 
uj 25 = 2a£n 2 , (55) 
Wj = 0, i = 26,...30. 

Thus, condition §§§ cannot be satisfied and Dkij do not 
evolve according to parabolic equations either. 

Therefore, system 1 is not parabolic and thus theorem 
1 cannot be used to prove well-posedness. However, the 
eigenvalues of the Fourier transformed second-order part 
of the evolution equations constitute a diagnostic tool for 
well-posedness and in order for system 1 to be well-posed, 
the following conditions are necessary 

A >,9 > O,0< 0,£ < 0,e > 0. (56) 

Before a sufficient proof of well-posedness is given, note 
that in appendix [B] it is shown that if the Fourier trans- 
form P(ik) of the principal part of a quasi-linear spatial 
operator P(d x ) 

a) satisfies the Petrovskii condition [33|, i.e., for all k 
the eigenvalues u> of P(in) satisfy Re{u>) < c, where 

c e R, 

b) it has a complete set of eigenvectors and 

c) the eigenvectors can be chosen as analytic functions 
for k^O 

then equation dtU = P(d x )u admits a well-posed Cauchy 
problem. 

For general second order systems which are not 
parabolic, the first order terms have to be taken into 
consideration in the definition of the principal part [29| . 
Hence, one has to consider both second and first order 
terms in the principal part of system 1. 

For the family of gauges (| 19[) with zero shift vector, 
the Fourier transformed principal part of system 1 has 
28 non-zero eigenvalues, 22 of which have a non-zero real 
part and the remaining 6 are purely imaginary. Due to 
the Petrovskii condition, only the real part of the eigen- 
values is important for establishing well-posedness. It 
turns out that the Petrovski condition can be satisfied 
only when 

£ = 4>- (57) 



If ([57]) is true, the real parts of the 22 complex eigenvalues 
are equal to the 22 non-zero eigenvalues ([53"]) - ([55]) of the 
second-order part of system 1. Hence, if both ([56"]) and 
([57]) hold true, the Petrovskii condition is satisfied with 
c = because the real part of all eigenvalues is strictly 
non-positive. 

The Fourier transform of the principal part of system 
1 is found to have a complete set of eigenvectors, if one 
of the following conditions holds true 

6 y£ 2e, A ^ e (58) 

or 

9 ^ 2e, A = e = -2£ = -20 (59) 

Finally, a calculation of the eigenvectors of P(in) of 
system 1 shows that the eigenvectors can be chosen as 
analytic functions for k =^ 0. In particular, their compo- 
nents can be chosen as sums and products of polynomials 
of Ki and polynomials of k. 

Therefore, if ([56"]) and ([57]) are satisfied and one of 
equations ([58]) or ([59]) is true, then conditions a), b) and 
c) above are met and system 1 admits a well-posed initial 
value problem. This holds true, if a non-zero fixed shift 
vector is assumed. 

Note here, that the presence of the lapse function as 
a multiplicative factor in the added terms of system 1 
is not necessary for well-posedness and omitting it may 
prove useful. This is because the damping rate of the 
constraints is of the form exp(— a\K 2 t), and hence when 
the lapse function is small, the damping time-scale may 
become large. Examples where the lapse function attains 
small values are black hole spacetimes. It is well-known 
that a becomes small near the singularities. A way to 
overcome the difficulty of inefficient damping, when the 
lapse function is small, is to allow large values for A. How- 
ever, this could result in numerical complications (see sec- 
tion [VlT]). Therefore, it may eventually be better to re- 
move the lapse function from the damping terms because 
this results in a damping rate of the form exp(— \n 2 t). 

Another important subject is the structure of the evo- 
lution equations of system 1. In addition to the complex 
eigevalues, the principal part of system 1 has 8 eigenval- 
ues with zero real parts. These are a subset of the eigen- 
values of the principal part of the ADM equations and 
are: with multiplicity two, ±ian each with multiplicity 
two and ±ia^/2An. Note that, if a non-zero shift vec- 
tor is considered, the two zero eigenvalues become f3 l Ki. 
These 8 modes correspond (as is illustrated in section lVl]) 
to constraint-satisfying degrees of freedom, which are the 
four gauge modes and the four physical modes. In [2jJ it 
was shown that for gauge (fl"9]l the constraint satisfying 
modes are described by strongly hyperbolic equations. 
Thus system 1 has a set of hyperbolic modes and a set of 
parabolic modes and for this reason it may be classified 
as a mixed hyperbolic - second-order parabolic (MHSP) 
set of equations. However, notice that it does not really 
fulfill the requirements of definition 2. 
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The goal of maximazing the number of non-zero eigen- 
values has been achieved by system 1. All degrees of 
freedom now evolve according to a principal part. This, 
according to 17[ , may be an advantage for numerical ap- 
plications. 

In addition, something even more important has been 
achieved: The parabolic modes correspond to constraint 
violating perturbations, which are now getting damped. 
This can be more easily realized, if one considers a so- 
lution to the equations and apply small amplitude con- 
straint satisfying and small amplitude constraint violat- 
ing perturbations separately. If constraint satisfying per- 
turbations are imposed on the system, then an exact can- 
cellation of the added terms of system 1 occurs and the 
system reduces to that of ADM. A thorough study of the 
evolution of constraint satisfying perturbations has been 
carried out in [27[ and hence this is not the subject of 
this work. If however, one imposes constraint violating 
perturbations on system 1 then the parabolic terms are 
turned on. As a result, the constraint violating perturba- 
tions and in particular the short-wavelength ones, have 
no other choice, but to decay exponentially with time. 
This will also be realized by the structure of the evolu- 
tion equations of the constraints in the following section. 



B. Evolution of the Constraints with System 1 

By virtue of system 1 and equations (|A1[) - (|A6[) , it can 
be shown that the constraint evolution equations with 
system 1 are 

d Ckij = a\~f ab d a d b C klj - aC,C ktj + a(X - e)j ab d a C bki:i 
- afrijdkH + a\{d kl ab )d a C bl3 
+ (d k a)X 1 ab C blJ , 



(60) 

d C k ii 3 = aej ab d a d b C k iij - C,aC k uj 

+ 2(d [k a) [e^'daCu^ + CC m + SrdjW] 
+ 2a^d [klij )d l] H + 2ea{d [kl ab )C u]l3 , 

(61) 

d H ~ -2^aj ab d a d b n - ae 7 sl l kj r l dsd a C kH3 , (62) 



d Mi ~ ^aB^dadbMi -~(9 + 4<j>)a~f ab did a M b , 

(63) 

where the exact evolution equations of C k ij and C k uj have 
been presented, but only the principal parts of the evolu- 
tion equation of Ti and Mi have been given ("~" means 
equal to the principal part). 

The eigenvalues of the Fourier transform of the prin- 
cipal part of system (j6"0|) - ([6"3")) are all non-zero and they 
equal the non-zero eigenvalues of the second-order part 
of system 1 ([53" ]) -([55 ]) . Thus, equations ([50 ]) -([55 ]) are 
second-order parabolic, when condition (|56|) is satisfied 
and hence they admit a well-posed Cauchy problem. 



The parabolic structure of the evolution equations of 
the constraints implies that the constraints are damped 
with time. In particular, the shortest the wavelength of 
perturbations which violate the constraints, the faster 
those perturbations are damped because the damping 
rate is of the form exp(— a\n 2 t). The new evolution 
equations of the constraints suggest that the surface of 
constraints is a local attractor. 

Another important property of system 1 is that the 
constraint violatiions are smoothed with time. This is 
due to the wavelength-dependent nature of the damping 
timescale. 

Furthermore, the structure and hence the smoothing 
and damping properties of equations (j60|) - (f63"|) is inde- 
pendent of the gauge condition or the spacetime solu- 
tion, since no specific gauge or spacetime geometry was 
assumed in the derivation of the new constraint propa- 
gation equations. 

Finally, equations (|rJ0|) . (l6"T]) and (|6"3"|) are completely 
decoupled to highest-order from the rest of the system 
and are manifestly second-order parabolic, when (|56p is 
satisfied. However, (|62l) is coupled to second-order to 

m- 

At this point, one may wonder if the aforementioned 
discussion on damping and well-posedness is relevant 
when equation (|34|) is enforced at the RHS of (|62|). This 
is so, because this enforcement turns (|62|) into third- 
order. Even in this case, it is straightforward to show that 
the eigenvalues remain the same and that the constraint 
propagation system remains well-posed. The proof of 
well-posedness of the constraint evolution system, with 
(IMl) enforced at the RHS of (|6"2")l , is exactly the same 
as the proof of well-posedness for parabolic systems pre- 
sented in 29]. For this reason it is not repeated here. 
Therefore, the properties of the evolution equations of 
the constraints are independent of whether (fM)) is en- 
forced or not. This is to be anticipated because all forms 
the constraint propagation equations (|60p - (|63p may ob- 
tain by use of (|34p are equivalent with each other and 
physical conclusions do not depend on the form of the 
equations. 



V. PARABOLIC MODIFICATIONS TO THE 
KST FORMULATION 

In this section system 2 is introduced. This system is 
based on the KST formulation as presented in 0. Thus, 
before going into the details of system 2, first the KST 
formulation is reviewed and the original gauge choice of 
the KST system is extended to the family (fT9")l . 



A. The KST formulation 

The KST modification to the first-order ADM formu- 
lation, as originally introduced, is a parameterized family 
of strongly hyperbolic formulations. 
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The KST modification is 

d t K l3 ={ADM) + pcrryH + iia~/ ab C a(ij)b , 
d t D kij ={ADM) + Tjorfk^Mj) + xajijM k , 



(64) 



In addition to these modifications to the RHS of the 
ADM system, Kidder, Scheel and Teukolsky used a den- 
sitized lapse (|2"2")) and a fixed shift vector /3\ 

In this work the lapse condition is extended to the 
algebraic family (fT9|) . The resulting KST system is es- 
sentially the same as the original one, but with one dif- 
ference: a must be replaced by the quantity A defined in 
equation (|20p . For illustration, the principal part of the 
second-order covariant derivative of the lapse is 



1 



ViVja ~ -aA^idjDm + d t D jkl ). 



(65) 



The RHS of (|65|) is exactly the same as the RHS it has 
in the original KST system |2j, but with a replaced by 
A. 

The 30 characteristic speeds of the principal part of the 
KST formulation in the generalized gauge (p~9|) are equal 
to {0, ±a, ±cia, ±c2a, ±C3«}. Eighteen of the speeds are 
0. Speeds ci, C2, C3 are different than those of the original 
KST formulation and are given by 



2A, 

i(r? - 4r)A -2x- 12A X - 677V), 

i (2 + 8p - n - 4pr? + 2 X + 8px ~ 2r^), 



(66) 



where the squares of the speeds is shown. For the sys- 
tem to be well-posed the RHS of equation (|66|) must be 
positive. 

The analysis of well-posedness of the KST system in 
the extended gauge is the same as that of the original 
one. Hence, if the RHS of equation (|66[) is positive, the 
new evolution system is strongly hyperbolic unless one of 
the following conditions is satisfied 



Ci = 0, i =1, 3 
ci = c 3 7^1, 
ci = c 3 = 1 t^c 2 . 



(67) 



Ensuring globally that conditions (fBT)) will not be sat- 
isfied is not as straightforward anymore due to the spa- 
tial dependence of A. However, there is a way around 
this difficulty. Notice that if speeds c 2 , C3 are set to 
unity, the strong hyperbolicity conditions are met. This 
can be done irrespectively of whether densitized lapse or 
"1 + log" slicing is chosen, by requiring that the following 
conditions be satisfied 



r) + 3% = 0, r\ — 2\ — drjxjj =8, 
2 - 7} + 8p - A-qp + 2 X + 8px ~ 2^ =2- 
Solving these equations parametrically one obtains 
5x + 8 1 



(68) 



18 X 



X + {0,-8/5} 



(69) 



V 



X 



(70) 



or 

6 5 

5' ^ = "6' - 5' 

The choice x = —8/5 in (|69|) renders the system weakly 
hyperbolic and hence ill-posed. The solutions above co- 
incide for the following values 



V 



6 



x 



(71) 



The generalization of the gauge constitutes the first 
extension of the original KST formulation. 



B. Evolution of the constraints of the KST 
modification to ADM 

By virtue of equations (|Alj) - (IA6j) , one can show that 
the evolution of the constraints with the KST formulation 
becomes 



d Cki 



d Q H ~ a( V -2 X - 2)y mn d m M n , 



(72) 

0. 

(73) 
(74) 



d Q Mi 



- 5 a(l 

1 , 
--a(l 



Ap)d i H + -a(l 



2i,) 1 ab 1 ls 3 l C ai 



2^)7 



ab Y k d k c s 



-a(l + 2A) 1 ab l sk d k C siabl 



(75) 

where the exact evolution equations of the C k ij , C k uj con- 
straints are presented, but only the principal part of the 
evolution of TL 1 Mi is given. The apparent discrepancy in 
these equations with the equations presented in [2| , arises 
because of the difference of the definitions of TC and C k uj ■ 
In particular, TL — 2Ti and C k uj — 2Ckiij, where 7i and 
CkUj denote the definition of those constraints in (2j. 

The characteristic speeds of the Fourier trans- 
form of the principal part of equations (|72 p - (|75|) are 
{0, ±C2a, ±C3a}, i.e., they are a subset of the characteris- 
tic speeds of the evolution equations, as expected. It was 
shown in Q that if the evolution equations of the dynam- 
ical variables are strongly hyperbolic, then the constraint 
propagation equations are also strongly hyperbolic. 

Although this property is really desirable, the KST 
formulation has a serious disadvantage. Equation (|72p . 
shows that Ckij evolve according to lower order terms. 
This implies that "static" modes correspond to the prin- 
cipal of l|72p . These low order terms of (fT2"|) are functions 
of the momentum constraints. This in turn means that 
any violations of the momentum constraints will directly 
lead to violations of C k ij ■ This is not desirable and there- 
fore damping of Ckij should be introduced in the KST 
system. 
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C. Mixed Hyperbolic- Parabolic KST Formulation 

In this section the generalized KST system of sec- 
tion lV Al is modified so that the resulting formulation has 
the desirable properties of strong damping and smooth- 
ing of Ckij ■ This is achieved by the same parabolization 
strategy that was applied to the scalar wave equation of 
section [II] In particular, the following modification is 
suggested 

d oll] =(KST) + \ ai ab d a c bl3 , 

(76) 

d a D kij =(KST) + (aC klJ , 

where (KST) stands for the RHS of t he ext ended KST 
formulation as that is given in section IV Al Equations 
([76]) together with the evolution equations of -fQj from 
equations (|64|) are called system 2 throughout this work. 

For A > the evolution equation of the three-metric 
is a parabolic equation (see section IV C lj) . Notice also 
the term QaCkij which has been added to the evolution of 
Dkij ■ This is important to couple the Ckij constraint vi- 
olating modes which arise from the evolution equations 
of the three-metric, to those which arise from the evo- 
lution equations of the Dkij variables. It is straightfor- 
ward to check that if A = 0, then there still are eighteen 
zero eigenvalues in the principal part of the formulation, 
whereas if £ = then there are twelve zero eigenvalues. 
As in equation ([5]), if both A > and £ > 0, strong damp- 
ing of Ckij occurs. The A term takes care of the short 
wavelength perturbations, whereas the £ term takes care 
of the long wavelength ones. 

Note that the proposed modification can be easily in- 
tegrated with the system of redefined variables, which 
was introduced in [2j. This is so, because in that sys- 
tem only Kij and Dkij were transformed. This may be 
advantageous, because the redefined system has shown 
better stability properties in numerical simulations than 
those of the KST formulation presented in section IV Al 

1. Well-posedness of system 2 

To prove well-posedness for system 2, note that the 
equations can symbolically be written as 

d olij = (KST) + aXj^dadbjij - a\j ab d b D aij , (77) 
d Kij = (KST), (78) 

d D kij = (KST) + (a(d klij - Dkij). (79) 

It can be readily realized that equations (|7"7| - (fT9")) have 
obtained the form of a MHSP system (|10[) with u = Yy 
and v — {K.^, Dkij}- The coupling terms are i?i2 = 
—aXj ab d a D b ij and R^i = Qadk^ij which are first order. 
The uncoupled (|77p is a second-order parabolic equation 



provided that A > and the uncoupled system of ([78)) 
and (|T9"|) is strongly hyperbolic. Therefore, the total set 
of equations is mixed second-order parabolic - strongly 
hyperbolic, bearing thus resemblance to the compressible 
Navier-Stokes equations. According to theorem 2, this 
system admits a well-posed Cauchy problem. 

The value of £ does not play any role in the establish- 
ment well-posedness, but it is really important for the 
stability of numerical simulations. If £ < 0, then at some 
point exponential explosion of the error of Ckij will take 
place. Thus, as it was explained in section UT1 C > is 
required for long-term and stable numerical simulations. 

Finally, note is that the lapse function multiplying the 
A and £ terms in equation (|76[) may be dropped and the 
system will still be well-posed. As is discussed in section 
IVIII this may be advantageous in spacetimes, where black 
holes are involved. 



D. Evolution equations of the constraints of the 
parabolized KST formulation 

The modifications introduced to KST in equation (|76[) . 
result into modifications of the evolution equations of the 
constraints. These modifications can easily be derived 
by virtue of equations (|A1[) - (|A6[) and the resulting con- 
straint propagation system is 

doCkij = a\~i ab d a d b Cki 3 - a(,C k i 3 

+ aX^daC^ + a\(dkl ab )d a C bl] - (80) 
- airYmMj) - axiijMk- 

d t H ~ (...) + a(C {1) H aUl + \r b {i) H l ')d a C blJ 

dtCkiij — (■ ■ ■), 
d t Mi ~ (. . .), 

(81) 

where (. . .) stands for the RHS of equations (|75|) - (|75|) and 
where ( 1 ')'H abl i and ( 4 )7i y are given in appendix lAl In 
the equations above, the exact RHS has been written for 
the evolution of Ckij and only the principal parts for the 
remaining constraint variables. 

Equation (j50)) resembles the form of ([S]) and the evolu- 
tion of Ckij is not left at the mercy of the low order terms 
as in the KST formulation because the principal part is 
not zero anymore. 

Equations (|5D|) and (fSTj) have obtained the form of 
MHSP just like the evolution equations of system 2. In 
particular, equation (|80p is parabolic, the uncoupled sys- 
tem (|8ip is strongly hyperbolic and the coupling terms 
are of first order. Therefore, the constraint propaga- 
tion equations with system 2 admit a well-posed initial 
value problem. In addition, due to the parabolic nature 
of equation J80|) . the Ckij constraints are damped and 
smoothed with time. 
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VI. STABILITY OF FLAT SPACE 

In this section the stability of fiat space is studied in 
conjunction with system 1, the first-order ADM and the 
KST formulations. 

System 2 will not be considered any further in this 
work so that the theoretical predictions of this section are 
made more relevant to the numerical simulations shown 
in [3l[ . A flat space stability analysis along with numer- 
ical simulations with system 2 will be considered else- 
where. 

The stability analysis is carried out by linearizing 
the evolution equations and the constraints about flat 
space. In particular small perturbations are applied 
on all dynamical variables as follows 7^ = Sij + 5"/ij , 
Kij = + SKij, Dkij — + SDkij- The linearized evolu- 
tion equations show how these perturbations evolve with 
time, while the linearized constraint equations serve as 
a tool for the classification of the individual modes in 
constraint satisfying and constraint violating. Both the 
evolution and the constraint equations are studied in the 
Fourier space in order to reduce the analysis to an eigen- 
value problem. 

In what follows to make the comparison among the 
three formulations appropriate the same gauge condi- 
tion is employed. This gauge condition has zero shift 
vector and lapse function of the form a — a (7) with 
A = d In a/d In 7 > 0. When linearization is carried out 
about flat space, "1+log" slicing cannot be distinguished 
from densitized lapse, because A = 1 for both "1+log" 
slicing and a densitized lapse with <j = 1. 

In addition to a specific gauge condition, only pertur- 
bations along the a^-direction will be considered. The 
results for any direction are qualitatively the same, but 
because they are rather complicated, here only a one di- 
mensional analysis is presented. 

The linearized evolution equations are 

d t Sjij = -2SKij + Xd x d x Sjij - Xd x D Uj , (82) 



dtSKij = -[d x SD ja + 6ijd x SD Hk + d x 5D ljl 

+ 5 u d x 5D kjk - (1 + 2A)S u d x 8D jkk 

- (1 + 2A)S lj d x 6D ikk - 2d x SD Uj ] 

+ ip(d x D {ij}1 - 5 1{l d x D aj)a ) 

- <><),, <(>,i),<)l\ 22 + d x d x 5K 33 ) 

+ -^{hjd x d x 5Ku + SudxdxSKtj) 

- 8Su5ij(d x d x SKn + d x d x SK 2 2 + d x d x SK 33 ), 

(83) 



dtSDu-j = - 2d x SK lj + T]{8i^d x Kj)i - Sm6j)id x K mm ) 

- x5i 3 {d x K 2 2 + d x K 33 ) 

+ (,8ij(d x d x 5D m i m — d x d x 5Di mm ), 
d t SD 2 ij = rj^^dxKj^ - 5 2 (i5 j ) 1 d x K mm ) 

+ X$i]dxKi2, +ed x d x 5D 2i]1 
d t SD 3ij = rjlyS^idxKj^ - 5 3 (i5j)id x K mm ) 

+ X$i 3 d x K 13 + ed x d x SD 3ij . 

(84) 

Equations (|5!?|) - (f8"4")) contain all three formulations 
(namely ADM, KST, system 1) for brevity, because they 
all have the ADM one as the underlying formulation. 
Just by choosing the free parameters one can go from 
one formulation to another. So, the ADM equations 
are obtained by setting all free parameters to zero, the 
KST equations are extracted from (fB"2")) - (f8"4")) by setting 
{A = e = (j> = 9 = £ = 0} and finally system 1 is ob- 
tained from (|82|) - (|84|) . by setting p = i/; = x = 'r] = 0. 
Note that one could study the system of equations l[8"5|)- 
([51| as a whole, but this analysis is complicated and will 
be the subject of a future work. 

The linearized constraint equations are 



oCuj = 


dxS^ij - dDuj = 0, 




oCuj = 


- SD 2ij = 0, 




5C 3 ij = 


- 5D 3ij = 0, 




5H = 


d x SD 2 i2 + d x 5D 313 - 


- d x 5D 122 - d x 5D 133 = 0, 


SMi = 


- d x 6K 2 2 - d x SK 33 


= 0, 


SM 2 = 


d x SK 12 = 0, 




SM 3 EE 


d x 5K l3 = 0. 





(85) 

The next step is to Fourier transform equations (|52"|) - 
(|85p by requiring a solution of the form 



5K l0 =A^e 4(K;r - wt) , (86) 
8D kij =D kij e^ KX - ut \ 

The effect of equation is the same as replacing a 
partial spatial derivative, d x , with ik and a partial time- 
derivative, dt, with iuj. The Fourier transform of equa- 
tions ([52T) - ([55]) is presented in appendix [Cl 

To study flat space stability in conjunction with the 
ADM, KST and parabolized-ADM (system 1) equations, 
all one has to do is solve the eigenvalue problem which 
results by the Fourier transform of (f8"2"l) and which can 
be written in the form 

Mv = ujv (87) 
where the transpose of the column vector v is 

v T = {711,712,713,722,723,733,^11,^12,^13,^22,-^23, 

-K33, -Dm, -D112, -D113, -D122, -D123, D133, -D211, -D212, -D213, 
D222, D223, -D233, -D311, -D312, -D313, -D322, -D323, -D333} and 
M is the matrix whose eigenvalues and eigenvectors one 
has to look for. 
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Number of modes 


Speed Vi 


Violates C^ij 1 - 


Violates W 


Violates M{! 


2 





X 


X 


X 


2 


±V2A 


X 


X 


X 


4 


±1 


X 


X 


X 


2 


±1 


X 


/ 


/ 


16 





/ 


X 


X 



TABLE I: Classification of the ADM flat space modes into constraint satisfying and constraint violating. The first column gives 
the number of modes. The second column gives the characteristic speed, which corresponds to the modes of the first column. 
The third, fourth and fifth columns have either a "x" or a "/" as values. A "x" means that the constraint is not violated by 
the specific modes, whereas a checkmark "/" implies that at least one of the components of the constraint is violated by the 
specific modes. 



One can check whether these modes violate or satisfy 
the constraints, via pure inspection of equations (|C4[) . If 
the components of a given eigenvector of M are such that 
it makes all constraints vanish, the eigenvector is said to 
be a constraint-satisfying mode. Otherwise, it is called a 
constraint-violating mode. 



The 30 x 30 matrix M and its eigenvectors will not be 
presented in this work for brevity, but they are available 
by the author upon request. What will be presented 
are its eigenvalues oji or characteristic speeds Vi — u>i/K, 
where appropriate. 



A. Flat space stability with the first-order ADM 
formulation 

When the ADM equations are employed, M has ten 
non-zero eigenvalues and twenty-six eigenvectors, if A ^ 
1/2, whereas it has twenty-four eigenvectors, if A = 1/2. 
A classification of these modes into constraint satisfying 
and constraint violating is presented in table [U 

No matter what the value of A is, there always are eight 
constraint satisfying modes (first three rows of table fl} . 
In addition to these eight modes, there always are two 
constraint violating wave modes (fifth row of table []}. 

Note here that the results shown in table U are valid 
only for A ^ 1/2. The case A = 1/2 has to be consid- 
ered separately. The results with A = 1/2, arc qualita- 
tively the same, but as already mentioned matrix M has 
twenty-four eigenvectors. This 4-fold or 6-fold degener- 
acy is disadvantageous for the ADM formulation, because 
these degenerate subspaces lead to linear growth of the 
constraint violating modes. 

To see this one has to understand the structure of the 
4-dimensional degenerate subspace (A ^ 1/2). A partic- 
ularly useful tool for doing this is a Jordan decomposition 
(JD) of matrix M. If A — 1/2, one can perform the exact 
same procedure to study the six-dimensional degenerate 
subspace. 

The JD is a similarity transformation M = XJX^ 1 , 
where X is a non-singular transformation matrix and J 
is called the Jordan matrix, which has the eigenvalues 
of M on its diagonal. If J is diagonal then M has a 
complete set of eigenvectors, if however J is not diago- 
nal, then the off diagonal elements (which can be placed 



either above or below the main diagonal) are equal to 
unity. These off-diagonal elements can give the informa- 
tion of the combination of dynamical variables, which are 
responsible for the fact that matrix M does not have a 
complete set of eigenvectors. This can be seen, if the 
eigenvalue problem ([57)) is written as follows 

J(X" 1 v) = u){X~ l v). (88) 

Equation (|88[) is an eigenvalue problem for matrix J. 
One can find the Jordan matrix J and the non-singular 
transformation matrix X by standard methods and hence 
one can also find which components of the transformed 
vector X~~ 1 v correspond to the off-diagonal elements of 
J. 

By carrying out the JD in the ADM case one finds 
that the four zero modes, which are thought to be "non- 
evolving" and which belong to the four-dimensional de- 
generate subspace, arc described by the following equa- 
tions 

ai2 = -2K 12 , (89) 
d t K 12 = --Ad x D 211 -^A8 x D 222 

+ 2A)d x D 233 + ^d x D 323l (90) 
9 t7l3 = -2K 13 , (91) 
d t K 13 = ^d x D 223 - ^Ad x D 311 

~(l + 2A)d x D 322 -^Ad x D 333 . (92) 

Plugging in the components of the zero-speed eigen- 
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Number of modes 


Speed Hi 


Violates Ckij 1 - 


Violates W 


Violates Mil 


2 





X 


X 


X 


2 


±72A 


X 


X 


X 


4 


±1 


X 


X 


X 


4 


±1 


/ 


X 


/ 


2 


±1 


/ 


/ 


/ 


14 





/ 


X 


X 


2 





/ 


/ 


X 



TABLE II: Classification of the KST flat space modes into constraint satisfying and constraint violating. The first column gives 
the number of modes. The second column gives the characteristic speed, which corresponds to the modes of the first column. 
The third, fourth and fifth columns have either a "x" or a "/" as values. A "x" means that the constraint is not violated by 
the specific modes, whereas a checkmark "/" implies that at least one of the components of the constraint is violated by the 
specific modes.. 



vectors shows that the RHS of (JSD]) and ^ is 0. Hence, 
components 712, if 13, 7i3, if 13 should not evolve with 
time. This is true only if the RHS is exactly zero. How- 
ever, small numerical perturbations can result in an non- 
zero RHS. As a result, if 12 and if 13 will begin to grow 
linearly in time. Non-zero values of if 12 and if 13, in 
turn lead to violations of .M2 and M3 respectively. This 
linear growth in the non-linear regime, where the pertur- 



bations are strong, is likely to obtain explosive behavior 
and eventually terminate simulations. This is precisely 
the root of the ill-posed nature of the ADM formulation. 

So far, even the analysis of flat space-time shows that 
one has to seek for other formulations, which at least a 
complete set of eigenvectors. This is exactly what was 
achieved by the KST formulation. 



B. Flat space stability with the KST formulation 

If the KST formulation is employed, matrix M of equa- 
tion ([87)1 has twelve non-zero eigenvalues and a complete 
set of eigenvectors, irrespectively of the value of A. Nev- 
ertheless, it still has sixteen zero eigenvalues which cor- 
respond to constraint violating solutions. 

A classification of the KST modes into constraint sat- 
isfying and constraint violating is given in table HH Note, 
that the results presented in table [Til are for the values of 
parameters given by equation (|69p . 

A few remarks are worth making at this point: a) 
The components of the constraint satisfying modes of the 
KST formulation (first 3 rows of tabic ITlj) are exactly the 
same as those of the ADM formulation independently of 
whether A = 1/2 or A 7^ 1/2, as expected. 

b) The KST formulation has cured the pathology of the 
incompleteness of eigenvectors of the ADM formulation. 
Note further that if x = — 2/5, the number of modes that 
violate the Hamiltonian constraint are minimized. It is 
not known yet, whether this is true for arbitrary space- 
times. For this, a detailed study is required, which is 
outside the scope of this work. However, it is interest- 
ing to notice the coincidence between this property and 
equations ([7T]) . 

c) Despite the desired property of strong hyperbolicity, 
the KST formulation has eighteen static modes sixteen 



of which violate the Ckij constraints. According to [l7j |. 
the existence of static modes can lead to exponential ex- 
plosion caused by low-order term perturbations. 

The disadvantages of the KST formulation were 
demonstrated in [311 ] , where it was shown that numerical 
simulations with the KST system in conjunction with fi- 
nite difference methods are unexpectedly terminated in 
certain cases. In particular, the KST simulations of gauge 
waves and expanding Gowdy spacetimes were terminated 
even earlier than the ADM simulations. In addition, 
linear growth of Ckij was observed in the evolution of 
initially small random noise. The unexpected behavior 
of the KST formulation was found to be caused by the 
growth of "static modes" , which were excited by the un- 
avoidable small violations of the momentum constraints. 



C. Flat space stability with system 1 

The flat space stability analysis with system 1, shows 
that system 1 is the most stable of the 3 formulations. 
For the choice of parameters that make system 1 well- 
posed, i.e. when both conditions (156|) and ([57]) are satis- 
fied, and one of equations (T58[) or ([59[) holds true, M has 
twenty-eight non-zero eigenvalues and a complete set of 
eigenvectors. 

The results presented in this section are specifically for 
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Number of modes 


Eigenvalue uii 


Violates Cm-,? 


Violates W? 


Violates Mi? 


2 





X 


X 


X 


2 


±« 


X 


X 


X 


4 


±K 


X 


X 


X 


2 


— i^£K 2 


X 


X 


/ 


10 


—ieK 2 


/ 


X 


/ 


8 


—ieK 2 


/ 


X 


X 


2 


—ien 2 ± k 


/ 


/ 


/ 



TABLE III: Classification of the flat space modes of system 1 into constraint satisfying and constraint violating. The first 
column gives the number modes. The second column gives the eigenvalues, which correspond to the modes of the first column. 
If the eigenvalue is purely real, then the corresponding mode propagates as a wave, whereas if the imaginary part of the 
eigenvalue is negative, the corresponding mode is getting damped. The third, fourth and fifth columns have either a "x" or a 
"/" as values. A "x" means that the constraint is not violated by the specific modes, whereas a checkmark "/" implies that 
at least one of the components of the constraint is violated by the specific modes. 



parameters which satisfy (fS"5|) and £ = 0. The analysis 
with parameters satisfying (|59[) and £ = is qualitatively 
the same. The reason why £ is set to zero, is because the 
results with non-zero £ are rather complicated. The ef- 
fect of the £ term has already been discussed both in sec- 
tion [IT] and in 0] and that is the wavelength-independent 
damping of all modes which violate Ckij ■ 

In addition to the choice of parameters described 
above, a harmonic lapse function (A = 1/2) has been 
used here for simplicity. The results remain qualitatively 
the same for any value of A. A classification of the sys- 
tem 1 flat space modes into constraint satisfying and con- 
straint violating follows in table 11111 Notice that in this 
case the eigenvalues have been given, instead of the char- 
acteristic speeds. 

A few remarks can be readily made: a) The eigenvalues 
of all constraint violating modes are non-zero and in par- 
ticular they have negative imaginary part, which implies 
exponential damping. The damping timescale is inversely 
proportional to k 2 , which means that the short wave- 
length perturbations are damped extremely efficiently. 

b) Another comment regards the nature of the sys- 
tem 1 modes. The first eight constraint-satisfying modes 
propagate as waves (the eigenvalues are purely real), ex- 
cept for two which are static. These modes are the same 
among all formulations, which are derived by addition 



of combinations of constraints to the RHS of the ADM 
formulation. 

All remaining modes are only damped, except for the 
two of the last row of table IIIIl These two modes repre- 
sent damped traveling waves. The origin of these modes 
can be understood, if one takes a closer look at the ADM 
modes. In table [1] it is seen that apart from the six 
traveling constraint satisfying modes, there are another 
two traveling wave modes which violate the constraints. 
The parabolic terms of system 1 are second order (see 
equations (|52"|) - ({51)) ) and therefore they cannot interact 
with the first order terms of the ADM formulation and 
lead to an exact cancellation of the two ADM travel- 
ing constraint-violating modes. However, because these 
ADM modes violate the constraints, system 1 damps 
them while they propagate. 

c) Finally, the zero eigenvalues (first row of table (JTTTJ) ) 
correspond to constraint satisfying modes and thus they 
do not pose any threat. Nevertheless, a fixed shift j3 l = 
constant assigns non-zero speed to these two modes. In 
particuarly they become equal to (3 l Ki, where Ki is the 
one- form which corresponds to the wavevector. 

All aforementioned predictions about the damping 
properties of system 1 have been numerically confirmed 
in [3l|. 



VII. ON THE CHOICE OF PARAMETERS 



This section discusses some ideas, which are related to 
the choice of parameters of system 1 and 2. It should be 
noted that this is a complicated discussion and only pre- 
liminary arguments will be made which are by no means 
complete. An analysis of this delicate subject has also 
been given in [3lj . 



Because the equations have been parabolized, the val- 
ues of the free parameters depend mainly on two factors: 
a) The effectiveness of the damping rate and b) the nu- 
merical scheme employed to integrate the evolution equa- 
tions. 

As already mentioned in sections |TT] and IIVB1 the 
damping rate is of the form exp(— a\n 2 t). Therefore, 
what one would desire is for the minimum damping rate 
to be efficient enough. The minimum damping rate is 
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exp(— a m i n \K^ nin t) , where a m i n is the minimum value of 
the lapse function and K m i n = 2-k/L with L the length 
of the computational domain. Therefore, if a m i n is very 
small then the damping rate is very small, unless the 
system 1 parameters are assigned large values. 

In principle, setting the damping parameters to very 
large values results in a very efficient drive towards the 
surface of constraints. Nevertheless, there are numerical 
issues involved with parabolic equations which have to be 
addressed. Before the validity of the following discussion 
is taken for granted a Von Neuman stability analysis [34] 
of the evolution equations is required. However, general 
statements can still be made even if such an analysis is 
not available. 

The numerical stability of an explicit numerical scheme 
for parabolic equations requires that the mesh size and 
damping parameter satisfy a Courant condition of the 
form XAt/ (Ax) 2 < c, where At is the time step, Aa; is the 
size of the mesh and c is a constant. For example, c = 1/2 
for a one dimensional linear diffusion equation with a four 
point explicit scheme with forward time difference [34|. 

The above Courant condition implies that the time 
step has to be really small, if A is large. It is this tradeoff 
between time marching versus strong damping, which is 
the limiting factor on the values of the damping param- 
eters with explicit numerical schemes. 

Similar constraints on the maximum time step exist 
when explicit schemes are chosen for the ADM evolu- 
tion equations. In this case, it is the speed of propaga- 
tion of the information, which sets the limit on the time 
step. For most applications, numerical stability with the 
ADM formulation requires a Courant condition of the 
form [i = At/ Ax < c, where fi is the usual Courant num- 
ber. The same is true for other hyperbolic systems. For 
example, the numerical integration of a one-dimensional 
wave equation with characteristic speed unity and a stan- 
dard explicit scheme [lij requires fx — At/ Ax < 1. 

Systems 1 and 2 contain both parabolic and hyper- 
bolic terms, and hence for numerical stability both the 
parabolic and the hyperbolic condition have to be satis- 
fied [34|. Combining the two conditions, one finds that 
for system 1 in conjunction with explicit schemes the 
overall numerical stability is of the form A/i < cAx. 
Thus, A should have a value, such that \i and Ax 
should be large enough in order to finish computations 
within reasonable time-scales, while maintaining satis- 
factory damping properties. This has been shown to be 
possible in [31] . 

Note, however, that the above limitation on the time- 
step is not present when the numerical integration is car- 
ried out with implicit schemes. For example, a Crank 
Nickolson numerical scheme for the one-dimensional dif- 
fusion equation 34] is always stable. An implicit scheme 
would allow one to set the damping parameters to arbi- 
trarily large values without having to worry about limita- 
tions on the grid structure, but at the cost of cumbersome 
programming. 

Taking all these facts into account, it is concluded that 



until more detailed numerical simulations become avail- 
able the safe choice for the parameters of systems 1 and 2 
are those which satisfy the conditions for well-posedness 
presented in sections llV Al and IV C II respectively. 



VIII. CONCLUSIONS 

Two new parameterized 3+1 formulations have been 
presented. These two formulations are referred to as sys- 
tem 1 and system 2. 

The well-posed system 1 is derived by adding combi- 
nations of the constraint equations and their derivatives 
to the RHS of the first-order ADM formulation. As a 
result the evolution equations of the dynamical variables 
have become second order in spatial derivatives and it has 
been shown that system 1 has obtained mixed hyperbolic 
- second-order parabolic nature. 

The evolution of the constraints with system 1 is gov- 
erned by 2nd-order parabolic PDEs irrespectively of the 
gauge condition employed. This in turn means that the 
constraints are damped and smoothed in time, and hence 
the surface of constraints is a local attractor. 

The well-posed system 2, is a parameterized parabo- 
lization of what is usually called the KST formulation. 
The original KST gauge condition was generalized and a 
certain combination of the terms of system 1 was added 
to the RHS of the KST formulation. System 2 is a man- 
ifestly mixed hyperbolic - second-order parabolic system 
of equations. 

The effect of the terms introduced in system 2 is the ef- 
ficient damping and smoothing of perturbations which vi- 
olate a subset of the constraints of the formulation, while 
the remaining constraints evolve according to strongly 
hyperbolic equations. 

As a first test for system 1, a study of the stability 
of flat space against small amplitude perturbations was 
carried out and contrasted to the stability of the ADM 
and the KST formulations. The results indicate that in 
this case system 1 is more stable than both the ADM 
and the KST formulations. In particular, it was shown 
that system 1 damps exponentially and smooths all fi- 
nite wavelength constraint- violating modes. Numerical 
support for this prediction has been given in [3l| . 

These two new formulations provide a new basis for 
obtaining well-posed formulations of GR which possess 
constraint damping and in turn a new basis for construct- 
ing numerical schemes for long-term and stable numerical 
simulations of the Einstein equations. 
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APPENDIX A: EVOLUTION OF THE 
CONSTRAINT EQUATIONS 

In this appendix the approach of [16| is followed and 
the evolution equations of the constraints of the first or- 
der ADM formulation (not of the standard ADM formu- 
lation as was done in [161 ]) are given in terms of the time 
derivatives of the dynamical variables instead of in terms 
of the spatial derivatives of the constraints and the con- 
straints themselves (cf. equations ([3"5[) - (HT1) ). This is very 
useful because the addition of constraints to the RHS of 



the ADM formulation has a direct impact on the evolu- 
tion equations of the constraints, which becomes straight- 
forward to calculate if the evolution equations of the con- 
straints are expressed in terms of the time derivatives of 
the dynamical variables. 

It is straightforward to check that by taking the time- 
derivative d of the constraint equations of the first order 
ADM formulation one obtains 

d Ckij = dk0o7ij) ~ {d D kij ), (Al) 
d CkUj = d k (d Diij) - di(d D kij ), (A2) 



d o n=^H ski ^d s (d D kij )+ WH k v0 o D kij )+ ^H^(d K l3 )+ WHV(d olij ), (A3) 

where 

(l)jjskij _ ^kj^si _ jij^ka 

Q 1 

y 'tl J = ~7 J 7 7 L 'mns - -jl 7 7 L>smn ~ *1 7 T^mns 

i s?' ,17171 7-j | -.mi-.nj „,lk t-i „M~, s i ~,km 7-> 

+ 7 7 J 7 iJ smn + -7 7 J 7 L>i mn -j 7 J 7 D srnh 

< 4 >ff« = 2K m K n i - 2R» - 2KK^ + 1 mn T\ nnl ks Ths ~ l %a J 3b r s ak T k sb 

- 7 "V b l sl C as u - l sl l m3 l ln l ab D amn D hsl + 7"7 nJ 7 m V fc A mfe C smfc , 

d M % = ^M kab d k {d K ab ) + ^Mi ab o K ab ) + ^M kab o Dkab) + (4) M t ab o 7a 6 ), (A5) 



where 

(i) M .fea6 = 7 ka §i b _ 7 ^> ; 

(2) M .a6 = 7 a„ r b. n _ 
(3) M .kab = ^ab K .k _ ^ka R b + Iftabfi k^ 
(4) M .a6 = V . K ab _ ^jaj^b^ + ^^"^.6 _ R an T b m . 

(A6) 

Note that quantities ^M kab , ( ^Mi ab , ^M^, ( ^Mi ab , 

(i) H ski^ m H Hj^ (3) H ij and (4) H ij herej should not 

be confused with the similar objects H mns , ( 1 )M rrasn , 
{2) M , msn l^ (3) M ™ n and (4) M .msn in equations ([38l) - (|4T]) . 

A final note concerns the WfP 5 object. By a closer 
look, one realizes that the constraint Ckiij is part of a 
term in the definition of ( 4 )if^ . This means that addi- 
tion of constraints to the RHS of the first-order ADM 
evolution equations, will result in evolution equations for 
the constraints, which will have quadratic dependence on 
the constraints. This resulted, because in the derivation 



of (|Aip - (|A6p the Ricci tensor has not been calculated 
through (|31)|) . If ([56]) is enforced, then the C k uj term 
vanishes from the expression of WH ij , but the expres- 
sions of the equations above have to change accordingly. 

APPENDIX B: WELL-POSEDNESS OF SYSTEMS 
WHICH SATISFY THE PETROVSKII 
CONDITION 

This appendix provides conditions for well-posedness 
of the initial value problem with quasi-linear systems of 
PDEs which satisfy the Petrovskii condition. 

Let 

d t u = P(u,di)u, i = l,...s (Bl) 

be a set of quasi-linear PDEs, where u is the column 
vector of n dynamical variables, P is a quasi-linear dif- 
ferential operator, di denotes a spatial derivative and s 
is the number of spatial dimensions. For well-posedness 
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P can be considered to coincide with its linearized and 
localized principal part [29|], which is denoted by P. 

Before the proposition of this appendix and its proof 
are presented it is important to introduce a definition 
and a theorem from [29(. 

a. Definition 3: An operator, P, is called semi- 
bounded if there exists a Hermitian operator such that 
(u,Pu)h + {Pu,u)h < 2c(u,u)h, for all sufficiently 
smooth functions u(x) and c € R is a constant. 

The quantity (u,w)h above is called an energy inner 
product and its exact definition can be found in [291 ] . 

b. Theorem 3: If the operator P is semi-bounded 
the Cauchy problem for (|Blj) is well-posed. 

Now the following proposition can be proven. 

c. Proposition 1: If the Fourier transform, P (in), of 
Pa) satisfies the Petrovskii condition, b) has a complete 
set of eigenvectors, and c) the eigenvectors can be chosen 
as analytic functions for all n 7^ 0, then equation (|Blj) 
admits a well-posed Cauchy problem. 

Proof: Proving well-posedness in the case n = is al- 
ways trivial because then P(in) — [29]. For n ^= 0, if 
P(in) satisfies the Petrovskii condition then the eigenval- 
ues u>i of P(in) satisfy 



Re{u>i) < c, 



(B2) 



where c 6 R and is independent of k. 

If P(ik) has a complete set of eigenvectors, then the 
matrix S(k)P{ik)S{k)~ x is diagonal with the eigenvalues 
of P(in) on the diagonal, provided that contains 
these eigenvectors as columns. 

Now, the hermitian matrix H(k) = S*S, where S* is 
the complex conjugate transpose of S, satisfies 

HP + P*H = S*SP + P*S*S, 

= s ( sps- 1 + (s*)- 1 p*s*)s, 



s* 



I LUl + LOl 



V : 



(B3) 



S. 



UJ r , 



< 2cH, 



where Ui denotes the complex conjugate of u>i . The defi- 
nition of inequalities between Hermitian matrices can be 
found in 

The energy inner product (u,w)h exists because of 
assumption c) [2^] of proposition 1. Then, the following 
energy estimate follows 

(u,Pu) H + (Pu,u) H = 

= (P(iK)u, H(k)u) + (u, H{n)P{in)u), 

u, (H(K)P(iK) + P{iK)*H(n))u^j , 

< 2c(u,u) H , 

(B4) 



where u is the Fourier transform of u. Parseval's relation, 
(u,w)h = (u,w)h, has been used in lines 1 and 3 and 
inequality (|B3[) has been used in going from line 2 to line 
3 above. 

Inequality (|B4[) proves that P is a semi-bounded oper- 
ator. Then according to theorem 3, (|B1[) admits a well- 
posed initial value problem. This completes the proof of 
proposition 1. 



APPENDIX C: FOURIER TRANSFORM OF 
LINEARIZED EVOLUTION AND CONSTRAINT 
EQUATIONS OF ADM, KST AND SYSTEM 1 

This appendix presents the Fourier transform of equa- 
tions ([82 ]) - ([85 ]) . 

The Fourier transformed evolution equations are 

ujjij — — 2i5Kij — i\K 2 ^ij + XnDuj, (CI) 

(jKij = - -n[Djn + SijDkik + £>iji+ 

+ hiDkjk - (1 + 2A)SiiD jkk - 

- (1 + 2A)5 Xj D ikk - 2D Uj ] + 

P^ij ^X-^mml ^lmm) 

- - £l(iAy>) 
+ i(l>K 2 6 l:j (K22+K33)- 

-i^0K 2 (6 lj K ll + 6 ll K lj )+ 

+ ie^SuS^iitu + k 22 + K33), 



(C2) 



+ K\8 il {K22 + K 33 ) - i^n Sij(D 

mlm 

u>D 2 ij = - K77((5 2 ( l -^j)i - 5 2 (i5j)iK mm ) 

- x$i]Ki2 - ien 2 D 2 i 3 , 

U)D 3ij = - Kr/^^Kj)! - S^iSj^Kmrn) 

- X$ijK 13 - itn 2 b ziv 

(C3) 

The Fourier transformed linearized constraint equations 
are 



Cuj = 




- D uj 


= 0, 


C 2 ij = 


D 2ij 


= 0, 






D 3ij 


= 0, 




it 


P > 212 


+ -D313 


- D 


Mi = 


K22- 


f ^"33 = 


0, 


M 2 = 


Kl2 - 


= 0, 




M 3 = 


K13 = 


= 0. 





(C4) 



Note that in [17[ it was shown that only three of the 
four constraints (Hamiltonian and Momentum) are inde- 
pendent, whereas in (|C4|1 all four are independent. The 
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apparent discrepancy arises, because the first-order ADM 
evolution equations are not reduced to equations for the 
three-metric only. This is exactly what was done in [TtJ , 
and hence the analysis presented in [13] was for a differ- 
ent set of equations than those of the first-order ADM 
formulation. 

The reason why this approach is not followed in this 
work is that the totality of the modes present in the 
first-order ADM formulation is larger than that of the 
second-order ADM formulation. The first-order ADM 
formulation has opened up an eighteen-dimensional sub- 



space, where more instabilities than the standard ADM 
formulation may arise. The totality of solutions of the 
first order ADM is greater than that of the standard 
second-order ADM and the two coincide, if and only if 
the Ckij constraints are satisfied [32| . The proposition 
in this work is to simulate the Einstein equations using 
variants of the first-order ADM formulation. Therefore, 
a study of the totality of the modes present in the first- 
order ADM formulation and its variants is required and 
the evolution equations will not be reduced to equations 
for the three metric only. 
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